; create a Surface density map from particle position (x,y,z) 
; to be used with sig2lens gravlens command

m=8.6e8
rhoc=2.7755e11

if (n_elements(read) eq 1) then begin
print, 'Data already loaded'
endif else begin

;readcol,'halo388599.out',x1,y1
;readcol, 'halo1402911.out',x2,y2
;readcol, 'halo15.out',x3,y3
;readcol, 'halo10.out',x4,y4,z4

readcol, 'halo10.out', x0, y0, z0

; r200 (Mpc/h)
; 0.863868     0.764250     0.914033     0.951089

read=1
endelse

ind=where(sqrt(x0^2+y0^2) lt 0.3 and sqrt(x0^2+y0^2+z0^2) lt 0.951)

; Gravlens needs positions in kpc/h (see gravlens sig2lens command)

x=x0[ind]*1e3
y=y0[ind]*1e3

print, "N=",n_elements(x)


xmin=min(x)
xmax=max(x)
ymin=min(y)
ymax=max(y)

nbins=500

xbinsize=(xmax-xmin)/float(nbins-1)
ybinsize=(ymax-ymin)/float(nbins-1)


;print, "Crop region: dx=", xmax-xmin, " dy=", ymax-ymin, " binsize=", xbinsize, "x", ybinsize

im = HIST_2D(x,y, MIN1=xmin, MAX1=xmax, MIN2=ymin, MAX2=ymax, BIN1=xbinsize, BIN2=ybinsize)

;use adaptative smoothing?

result=asmooth_fft(im,4.0,6.0,plot=0)

; write data to file, this format can be read from gravlens sig2lens command

openw, 1, 'sigma.txt'
printf, 1,  nbins,' ', nbins
for j=0,nbins-1 do for i=0,nbins-1 do printf, 1, result.asm_im[i,j]*m  
close,1

openw, 2, 'sigma.range'
printf, 2, xmin,' ', xmax
printf, 2, ymin ,' ', ymax
close,2

; plot results

!p.charthick=2
!p.color=0
!p.background=-1

loadct,3

set_plot, 'ps'

;window,0,retain=2

device,/encapsulated,filename='sigma_smooth.eps',/inches,xsize=5.3,scale_factor=1.0,font_size=14,/bold,/color,bits=16

tvscl,alog10(result.asm_im+0.1)

device,/close

device,/encapsulated,filename='sigma_smooth_surf.eps',/inches,xsize=5.3,scale_factor=1.0,font_size=14,/bold,/color,bits=16

shade_surf,result.asm_im

device,/close

device,/encapsulated,filename='sigma_cotour.eps',/inches,xsize=5.3,scale_factor=1.0,font_size=14,/bold,/color,bits=16
contour,result.asm_im,nlevels=15  
device,/close

end

